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TECHNIQUE FOR PREDICTING HIGH-FREQUENCY STABILITY CHARACTERISTICS 

OF GASEOUS-PROPELLANT COMBUSTORS 
by Richard J. Priem and Jefferson Y. S. Yang* 

Lewis Research Center 

SUMMARY 

A technique for predicting the stability characteristics of gaseous -propellant com- 
bustors is developed based on a model which assumes that the system is driven by cou- 
pling between the flow through the injector and the oscillating chamber pressure. The 
theoretical model uses a lumped parameter approach for the flow elements in the injec- 
tion system plus wave dynamics in the combustion chamber. Stability characteristics 
(frequency and decay or growth rates) were calculated for various combustor design and 
operating conditions to demonstrate the influence of various parameters on stability. 
These results show that the stability of a given combustor is determined by the oxidant to 
fuel mixture ratio and that design changes in the oxidizer side of the system have a much 
larger influence on stability than similar changes in the fuel system. 


INTRODUCTION 

Recent interest by NASA in using hydrogen -oxygen thrusters for the Space Shuttle 
attitude control propulsion system (ACPS) has resulted in an extensive technology pro- 
gram. In this program (ref. 1) the gas/gas feed system received the greatest amount of 
attention and technology effort. The gas/gas feed system offers the advantages of ver- 
satility, flexibility, and light weight and the ability to be developed into a reliable high 
performance, fully reusable system with excellent thruster pulsing performance (ref. 2). 
To achieve the desired reliability and reusability will require that the system be de- 
signed, tested, and proven to have the same '’dynamic" stability required for the Space 
Shuttle main engines (SSME) in the Space Shuttle Orbiter (ref. 3). 


a|c 

Summer Faculty Fellow at the Lewis Research Center in 1971 and 1972; Assistant 
Professor of Engineering, Pacific Lutheran University, Tacoma, Washington. 


For liquid/liquid and gas/liquid feed systems, as used in the SSME, several analyt- 
ical models (refs. 4 to 7) are available for predicting combustor stability characteristics. 
These models have been used extensively in engine development programs to ensure that 
preliminary designs had the desired stability. The object of the program reported herein 
is to provide an analytical model for gas/gas rockets to predict stability characteristics. 
The analytical model of reference 6 predicts that gaseous flow variations through the in- 
jector are responsible for many of the stability characteristics observed in gas/liquid in- 
jectors. Therefore, this model was used as the basis for an all gaseous -propellant 
system. 

The lumped -element model (ref. 6) for the dynamic flow characteristics of the gas- 
eous injector was used in this investigation. Since dynamic stability (ability to damp a 
high amplitude disturbance in a finite period of time) would be required of any engine us- 
ing gas/gas injection, the technique of solving for a neutral stability design point as used 
in references 5 to 6 was not considered adequate. Therefore, the analytical model was 
set up to calculate the frequency and decay rate (or growth rate if the engine is inherently 
unstable) for a specific combustor design and operating condition. This allows the de- 
signer to calculate the decay rate for his combustor to ensure that it meets the require- 
ments for a "dynamically" stable engine. 

After the analytical model was developed it was used to calculate the stability char- 
acteristics of a "standard" engine that might be used to meet the requirements for an 
ACPS thruster. Calculations were also performed for various perturbations of combus- 
tor design and operating parameters to demonstrate the usefulness of the model and the 
sensitivity of the stability characteristics to these parameters. To enable others to use 
the analytical model a complete listing of the computer program used to make the cal- 
culations, along with a description of the input and output and a sample calculation, is 
presented in the appendixes. 


SYMBOLS 

A q injector orifice area 

. A t nozzle throat area 

a chamber speed of sound 

‘BpBg eqs. 15(a) and (b) 

C characterizes mass "capacitance time" of injector dome (eq. 10(a)) 

injector orifice coefficient 
C * nozzle throat choked speed of sound 
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nozzle acoustic admittance, see eq. (16) 
gravitational constant 

characterizes flow "inductance" of injector duct (eq. 10(c)) 
n tn order Bessel function 
eq. (17) 
length 

longitudinal wave mode number, Z =0,1,2, etc. 
chamber Mach number 
molecular weight 

transverse wave mode numbers, m = 1. 84, 5.33, 8. 53, etc. for 
m = 3. 05, 6. 70, etc. for n = 2 

propellant burning response, see eq. (19) 

chamber flow response, see eq. (18) 

injector flow response, see eq. (9) 

oxidizer to fuel flow ratio 

pressure 

dome pressure for choked flow as defined in eq. (3) 

characterizes flow "conductance" of injector orifice (eq. 10(b)) 

chamber radius 

universal gas constant 

chamber radial direction 

complex frequency, a + iu> 

total temperature of propellant 

time 

dome volume of propellant 
chamber axial velocity 
mass flow rate of propellant 
mass flow rate through nozzle 
total mass flow rate 
chamber axial direction 


n= l; 



a oscillation growth rate or decay rate if negative 
y ratio of specific heats 

7j c * C * efficiency 
0 chamber tangential direction 

p density 

. r delay time 

co angular frequency 

co Q natural frequency 

Subscripts: 
d dome 

c chamber 

o orifice 

Superscripts: 

_ mean 

' perturbed, X - X 


THEORY 

The analytical model treats the rocket as a system consisting of the injector, com- 
bustion chamber, and a combustion region with the appropriate boundary and compat- 
ibility conditions at the interface of the various regions. The flow in each region re- 
sponds to an impressed acoustic pressure oscillation originating in the combustion cham- 
ber. These flow responses are influenced by the geometrical and gas dynamic conditions 
in the injector and combustion chamber. All dependent variables are written as the sum 
of a constant mean term and a small magnitude term that is harmonic in time. All gov- 
erning equations are thereby linearized by this small perturbation technique. 


I njector 

The flow response in the injector is based on the analysis of Feiler and Heidmann 
(ref. 6) with modifications to include compressible flow through the injector orifices. A 
schematic of the injector with its various elements is shown in figure 1. All the dimen- 
sions of the various elements in the injector are considered to be small compared to the 
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wavelength of the oscillations. This permits a lumped-parameter treatment of the var- 
ious elements in the injector. 


Propellant Supply Line and Dome 

A continuous flow of propellant to the supply dome is assumed. Since the flow is > 
compressible and isentropic, pressure oscillations will affect the instantaneous total 
mass in the dome so that the dome acts as a capacitor for the flow. Perturbing the dome 
total mass balance, we get 


Pd V d 

^d 



( 1 ) 


The dome mean pressure P d is determined by combining the mean pressure drop 
across the injector orifice with the mean chamber pressure: 


W = A C H P H *1 

° d d _ 1 r q T 


/rj \2/y /b \(y+D/yl 


( 2 ) 


and P d is solved using a curve -fitting technique and P Q = P c> 

For certain chamber pressures and propellant flow rates the flow may become 
choked. The dome pressure for choked flow is then given by 


P 


choke 



( 3 ) 


If the flow is actually choked for the given flow rate W, the dome pressure is given by 


P d = 


W 


( 4 ) 
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The dome gas density is described by the perfect gas equation 


P d = _ (5) 

d R T/M 
o w 

Injector Flow Duct and Orifice 

With a short length duct the flow can be assumed incompressible with uniform mean 
pressure within the duct. The perturbed momentum equation is used to obtain the pres- 
sure drop across the duct: 



( 6 ) 


To obtain the flow perturbation through the injector produced by the pressure perturba- 
tions on either side of the orifice the compressible orifice flow equation (2) is perturbed 
to obtain 
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Injector Flow Response 

The flow oscillation of the injector system is a function of the pressure oscillation in 
the combustion chamber. This function is expressed as an admittance as given by 
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( 8 ) 


[injector exit 

Substitution of equations (1), (6), and (7) into equation (8) gives 

Cs 


N. 


ln J o 

1 + CRs + CIs^ 


(9) 
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(10a) 


(10b) 


WL 


I = — 


P cV 


(10c) 


The real part of N in j indicates the degree to which the flow responds in phase with 
the impressed pressure oscillation. The imaginary part of N^. indicates the amount of 
flow oscillation out of phase with the impressed pressure oscillation. 


Combustor Chamber 

The mean flow in the combustion chamber is assumed to be uniform, one- 
dimensional, and inviscid. The gas in the combustion chamber is assumed to have the 
properties associated with the products of combustion for the mixture ratio being 
metered to the chamber. The gas properties (y, M^, and C ) as functions of mixture 
ratio were obtained from the tables for hydrogen -oxygen propellants in reference 8. The 
mean chamber flow conditions were then calculated from the following: 

Pressure: 



(11a) 


Speed of sound: 


a = yC’ 


v(y+l)/(y-l) 


( y + 1) 


(lib) 
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Density: 


Pc = 


yP c g 


(11c) 


Mach number: 


W, 


M = 


_ 

p a7rR 
H c c 


(lid) 


The natural frequency of the chamber was calculated assuming hard walls as follows: 


a */ 2 r R c Z> 

••xV" 


(12) 


The three -dimensional perturbed pressure and velocity fields in the chamber have 
previously been determined by Priem and Rice (ref. 9). With u> replaced by s/i, 


K = _ y 


[j n (mr)e W e st ] [s(e BlZ + Ke B2Z )+ M^e 5 * 2 + B 2 Ke B2Z )] (13) 
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(14) 


where 


B l = 


sM + V(sM) 2 + (1 - M 2 )(m 2 + s 2 ) 
1 - M 2 


(15a) 


B 2 = 


sM - V(sM) 2 + (1 - M 2 )(m 2 + s 2 ) 
1 - M 2 


(15b) 


and v„ is the perturbed axial velocity. A constant K is determined from the boundary 

z 

conditions as follows: 
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( 16 ) 


then 


K 


B t (l - M 2 + yG N M 2 ) + sM(yG N - 1) 


5 (l - M 2 + yG N M 2 ) 


(17) 


+ sM(yG N - 1) 


The flow perturbation response at the injector end of the chamber (z = -L ) can then 
be written as 


N c = 


W'P 


WP' 


"BjL -B 2 L 
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I + ysMle 
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(18) 


This flow perturbation response must be matched to the response produced by the 
injector -combustion process combination. 


Combustion Process 


The burning process which embodies the effects of propellant mixing and chemical 
reaction is assumed to be characterized by delay times r for the fuel and oxidizer. It 
is also assumed that the burning process occurs in a very thin region (relative to the 
length of the chamber) immediately downstream of the injector. 

The propellant burning response is assumed to be the sum of the oxidizer and 
fuel response functions. The individual oxidizer and fuel responses are weighted by their 
fractional mass flow rates with a delay time to obtain the following: 





WNinj * 5 


Ts\ 

/oxidizer 



(19) 


Matching the propellant burning response function with the chamber response function 
produces the following boundary condition for the chamber -injector interface: 

N b -N c =0 (20) 
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Equations (9) and (18) to (20) are solved together to obtain the chamber and injector flow 
responses along with the complex frequency. The imaginary part of the complex fre- 
quency describes the period of oscillation and the real part describes the damping rate. 
A positive value for the real part of the complex frequency means the system is sponta- 
neously unstable and the oscillations will grow in amplitude with time. 


Numerical Solution 

The flow chart for the computer program to calculate stability is shown in figure 2. 
The program listing, the program input formats, and a sample calculation are given in 
appendixes A to C. The solution procedure is to first determine the mean chamber con- 
ditions and to test if the injector flows are choked. A guessed value of the complex fre- 
quency s is then used to initiate an iterative scheme to converge on a consistent s 
which satisfies the compatibility conditions of equation (20). The iterative process re- 
duces the error between injector and chamber flow responses, which is the left side of 
equation (20), to zero by Taylor’s formula for two variables. 

For any given engine configuration and flow condition, there exist multiple solutions 
of the complex frequency s. These solutions can be determined from a table of flow re- 
sponse errors as functions of a range of s values. The lowest error values on the map 
will be at or near a solution. This point can be used as the assumed frequency to initiate 
the iterations or to check a solution. The quandrants on the complex plane in which the 
errors lie are also calculated as an additional check on the existence of a solution at the 
lowest point. A sample error map is tabulated in appendix C for frequencies between 
51 000 and 90 000 radians per second and growth rates of 750 to -650 reciprical seconds. 


ANALYTICAL RESULTS 

To examine the influence of injector and combustor parameters on stability, a base 
engine about which a parametric study could be performed was established. For the base 
engine it was assumed that the physical dimensions and mass flow rates would be repre- 
sentative of engines tested in the Space Shuttle attitude control propulsion system (ACPS) 
technology program. It was also assumed that the engine would have neutral stability (a 
disturbance would neither grow or decay) and that the oxidizer and fuel flow oscillations 
through the injector would be 180° out of phase with the chamber pressure oscillations. 
Furthermore, it was assumed that the ratio of the oxidizer to fuel flow oscillation was 
the same as the ratio of the oxygen to fuel flow rate. It was also assumed that the delay 
times of the fuel and oxidizer corresponded to a half period of the oscillations. 
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TABLE I. - STANDARD ENGINE 


(a) Combustion chamber 


Length, L c , ft 

0.5 

Radius, R c , ft 

0. 1575 

Nozzle throat area, A t , ft 2 

0.0185 

Nozzle throat choked speed of sound efficiency, t)q* 

1.0 

Nozzle admittance, G R 
Pressure wave numbers: 

0.9166 + 0.0 i 

For n = 0 

0 

For m = 1. 84 

1 

For 1=0 

0 

Standard growth rate , a, 1/sec 

0.0 

Standard angular frequency, w, Hz 

9950 

Standard chamber response, N c 

0.9174 + 0.0 i 


(b) Injector 



Oxidizer 

Fuel 

Q 

Dome volume, V^, ft 

0.0036574 

0.0024363 

Orifice length, L , ft 

0. 00027833 

0.0098092 

Orifice area, A q , ft 

0.0019514 

0.0030957 

Orifice coefficient, 

1.0 

1.0 

Mass flow rate, W, lb/sec 

2.658 

0.505 

Specific heat ratio, y 

1. 36 

1.41 

Temperature per molecular weight, 

16.875 

393.0 

T/M w , °R/lb 

Combustion delay time, r, sec 

0.000050255 

0. 000050255 


The physical and gas dynamical properties of the standard engine are shown in 
table I. Additional information on the standard engine can be found in the sample cal- 
culation in appendix C . 


Stability Characteristics of Standard Engine 

The stability of the standard engine operating at various fuel and oxidizer flow rates 
is shown in figure 3(a). The lines of constant growth rates (a/u>, fraction/cycle) orig- 
inate at the origin of the plot and correspond to constant values of O/F (oxidant flow 
rate/fuel flow rate). Examining the equations that are used in the solutions reveals that 
flow rates can be eliminated in the equations by dividing by total flow rate and converting 
the equations to O/F. Therefore, the engine O/F uniquely defines the stability of an 
engine, independent of flow rates, for fixed values of the other parameters. 
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Stability of the standard engine as a function of O/F is shown in figure 3(b). At 
very low and high mixture ratios the engine is stable. Between a mixture ratio of 0. 35 
and 5. 26 the engine is unstable. A growth rate of 0. 1 fraction per cycle is very unstable 
as an oscillation with an amplitude of 1 percent of chamber pressure would grow to an 
amplitude of 100 percent in 48 cycles or 5 milliseconds. 

At the very low O/F the fuel side of the injector is operating in the choked flow re- 
gime. Therefore, it does not respond to the pressure oscillations. The oxidizer side of 
the injector is operating at a very low injector pressure drop to chamber pressure ratio, 
which results in a large flow response, but the fraction of total propellant flow that is be- 
ing oscillated is small so the engine is stable. At the high O/F ratio the opposite is 
true with the oxidizer flow being choked and the fuel only being a small fraction of the 
total flow. In the intermediate O/F region both the fuel and oxidizer respond to a pres- 
sure oscillation with sufficient magnitude to make the engine unstable. 


Influence of Design Parameters on Stability 

To demonstrate the effect of the various design variables on stability, calculations 
to determine the growth or decay rate of the first transverse mode (n = 1, m = 1. 84, and 
l = 0) were made in which each design variable was individually increased and decreased 
about the value it had in the standard engine. The results are plotted in terms of growth 
rate as a function of the oxidant to fuel flow ratio (O/F) in figure 4. The range of each ■ 
variable was arbitrarily selected to represent a reasonable range over which the design 
might be changed. 

The influence of the chamber radius and nozzle throat area is shown in figures 4(a) 
and (b). Increasing either the chamber radius or throat area significantly improved sta- 
bility (curves are lower and less area under the curves). The reasons for the improved 
stability with these two variables are entirely different. The change in chamber radius 
changed the frequency of the oscillation with a resultant detuning of the system. A nozzle 
throat area increase decreased the chamber pressure and as a result increased the ratio 
of injector pressure drop to chamber pressure, thereby decoupling the injector flow from 
the chamber pressure oscillations. 

The influence of changes in oxidizer injector length, dome volume, and orifice area 
on stability is shown in figures 4(c), (e), and (g). The stability at high O/F ratios was 
not influenced by changing any of the oxidizer variables. This is because at high O/F's 
the oxidizer flow is choked and, therefore, does not respond to any pressure oscillations; 
therefore, the oxidizer variables are not important in stability under these conditions. 
Increasing the oxidizer injector length to 0. 03 feet made the engine stable over the entire 
mixture ratio. This is because the high inertial effect of a long orifice prevents the ox- 
ygen flow from responding to a pressure oscillation. Decreasing the oxidizer orifice 
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area had a similar effect. Decreasing the oxidizer dome volume improved stability by 
reducing the reservoir that stores the propellant when the flow oscillates. 

The influence of changes in the fuel variables are shown in figures 4(d), (f), and (h). 
The fuel variables did not influence stability at the low O/F ratios because in this re- 
gion the fuel flow is choked and does not respond to any pressure oscillations; thus, the 
fuel variables are not important under these conditions. At high O/F ratios the fuel 
side variables influenced stability in a manner similar to that described previously for 
the oxidizer side; however, the influence of the fuel variables on stability was much less 
than that observed with the oxidizer variables. 

Increasing the fuel and oxidizer flow properties (ratio of temperature to molecular 
weight) improved stability in a manner similar to those described previously for the in- 
jector (figs. 4(i) and (j)). Fuel properties had no influence on stability at the low O/F’s, 
and the oxidizer had no influence at high O/F's. Increasing the temperature over molec- 
ular weight ratio improved stability by increasing the injector pressure drop, thereby de- 
creasing the flow oscillations that could be produced by a given pressure oscillation. 
Again, changing the oxidant properties produced a larger influence on stability than did 
the fuel property changes. 

The influence of the oxidizer and fuel delay times on stability is shown in figures 4(k) 
and (1). Decreasing the oxidizer delay time to 3x10" 5 second (corresponding to burning 
in 0. 3 in. ) produced a very stable engine over the entire O/F region. Increasing the 
oxidizer delay time produced a stable operating region from an O/F of 0. 8 to 5. Again, 
a change in oxidizer delay time did not influence stability at the high O/F ratios. 
Changing the fuel delay time had a similar influence on stability as the oxidizer delay 
time but produced a much smaller change in stability. 

Looking at all the parameters together we see that the standard engine could be made 
very stable by increasing the chamber radius, throat area, oxidizer injector length, or 
reducing the oxidizer orifice area and delay time. Of these, the changes in throat area 
and orifice area influence the supply pressure and chamber pressure which might not be 
satisfactory from the overall systems point of view. Stability, therefore, would be eas- 
iest to obtain in this engine by increasing the oxidizer orifice length to 0. 05 feet (0. 6 in. ) 
and decreasing the delay time by a factor of 2 to 0.000025 second. 


SUMMARY OF RESULTS 

A technique for predicting the stability characteristics of gaseous propellant com- 
bustors has been developed based on a model which assumes that the system is driven by 
coupling between the flow through the injector and chamber pressure oscillations. The 
technique was used to calculate the influence of various combustor design and operating 
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parameters on stability. The results of these calculations may be summarized as 
follows: 

1. The stability of a given combustor was determined by the oxidant to fuel mixture 
ratio for any fuel or oxidant flow rate. 

2. Changing design parameters in the oxidizer side of the injector influences the sta- 
bility characteristics at low oxidant to fuel mixture ratios, while changes in the fuel sys- 
tem influenced the stability characteristics at high oxidant to fuel mixture ratios. 

3. Changes in the design of the oxidizer system had a much larger influence on sta- 
bility than a similar change in the fuel system. 

4. Changing the chamber radius or nozzle throat also had a large effect on the sta- 
bility of a combustor. 

5. To obtain maximum stability in a combustor with given propellant flow rates the 
combustor should have a large chamber radius and throat area. The injector oxidizer 
orifices should have a small total area, a long length, and should produce a very small 
delay time between when the oxidizer is injected and burned. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, June 25, 1973, 

502-04. 



APPENDIX A 


COMPUTER LISTINGS 


$ IBFTC MAIN 

REAL LL, MW , LOX , LF , IOX , IFU , MACH , NBR , NBI ,NCR , NCI 

COMPLEX CSROX , CSRF , ISOX , I SF , TAUSOX , TAUSF 

COMPLEX NC (3) ,SZ ,SZHZ ,GNOZ ,SS ,SSAVE j 

COMPLEX S (3) ,NB (3) ,NBOX,NBF ,ERR(3) , ESAVE1 , ESAVE2 

INTEGER WHICH, FRCUT,FICUT, REGION 

DIMENSION El (3) ,ER{ 3) 

DIMENSION F (200) ,T (15) 

DIMENSION Z1 (26) ,Z2 (26) ,IZ1 (26) ,IZ2 (26) ,Z3 (26) ,Z4 (26) ,IZ3 (26) 
COMMON/BF/F 

COMMON/ERM/FRMAX , FIMAX , FRMIN , FIMIN , FRCUT , FICUT , IOX , ROX , COX , 
*TAUOX , CFU , RFU , IFU , TAUF , IFLAGO , IFLAGF , WF , WOX 
COMMON/WH/ R , AT , TOMOX , TOMF , AOX , AF , LOX , LF , VOX , VF 
COMMON /SCLL/VZ , AMOR,AMV ,GC , W ,LL , AA , GNR, GNI 

COMMON /AMAIN/PI, GCON,GASC, PC, WTOT 
READ (5,150) WMAX , WMIN , GAMOX , CORFOX , AOX ,LOX 
READ (5,149) VOX , TAUOX , TOMOX , IP ,DELX 
READ (5,150) WFMAX ,WFMIN ,GAMF ,CORFF ,AF ,LF 
READ (5,149) VF , TAUF , TOMF , MAP 
READ (5,151) LL , R, AT ,CSTREF , GNR , GNI , SM , SL 
READ (5, 174) WOSTD , WFSTD ,DW, DWF 
IF (IP .EQ. 0) READ (5,175) FRINP ,FIINP , INPUTF 
IF ( IP . NE . 0 ) READ (5,174) (Z1 (I) , Z 2 (I) , Z3 (I) , Z4 (I) , IZ1 (I) , 

*IZ2 (I) , IZ 3 ( I ) ,1=1, IP) 

IF (MAP . NE . 0 ) READ (5,174) FRMAX , FRMIN , F IMAX , FIMIN , FRCUT , F ICUT 

IF (IP.NE.0) GO TO 1040 

IP=1 

Z 3 ( 1) =FRINP 
Z4 ( 1 ) =FIINP 
IZ1 (1)=0 
IZ3 (1) =1 

1040 DO 1000 MN=1 , IP 
IM=MN-1 

IF(MN.EQ.l) IM=1 
CALL RSTRP (MN , IZ1 (IM) ) 

1001 PMAX=Zl (MN) 

PMIN=Z2 (MN) 

FRINP=Z3 (MN) 

FIINP=Z4 (MN) 

WHICH=IZ1 (MN) 

NP=IZ2 (MN) 

INPUTF=IZ3 (MN) 

GNOZ=CMPLX (GNR, GNI) 

WRITE (6,152) 

WRITE (6,171) 

WRITE (6,153) 

WRITE (6,154) 

WRITE (6,155) WMAX , WMIN , GAMOX , CORFOX , AOX , LOX , VOX , TAUOX , TOMOX , DW , 
*WOSTD 
WRITE (6,156) 

WRITE (6, 154) 

WRITE (6,155) WFMAX , WFMIN , GAMF , CORFF , AF , LF , VF , TAUF , TOMF , DWF , WFSTD 
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WRITE (6, 157) 

WRITE (6,158) 

WRITE (6,155) LL,R,AT ,CSTREF ,SM, SL ,GNOZ , IP 
WRITE (6, 710) 

WRITE (6,155) DELX , PMAX , PMIN , WHICH , NP , MAP 
WRITE(6, 730) 

WRITE (6,155) FRINP ,FIINP , INPUTF 

FIINP=FIINP*2.*PI 

WRITE (6,172) 

DP=0. 

IF(NP.EQ.l) GO TO 401 
DP= (PMAX-PMIN) /FLOAT (NP-1) 

401 DO 500 JP=1 ,NP 
REGION=0 

IF (WHICH. EQ.O) GO TO 440 
CALL WHICHP (WHICH, JP , PMAX, DP) 

440 REGION=REGION+1 
WOX=WOSTD 

GO TO (400,410,410,400), REGION 
400 NO=IFIX (ABS (WOSTD-WMIN) /DW) +1 
GO TO 415 

410 NO=IFIX (ABS (WOSTD-WMAX) /DW) +1 
415 IF ( REGION. LE. 2) GO TO 3 

NO=NO-l 
WOX=WOSTD-DW 

IF (REGION. EQ. 3) WOX=WOSTD+DW 
3 DO 70 M=l,NO 

GOTO (420,430,420,430) , REGION 
420 NF=IFIX (ABS (WFSTD-WFMAX) /DWF) +1 
GO TO 250 

430 NF=IFIX (ABS (WFSTD-WFMAX) /DWF) +1 
250 WF=WFSTD 

IF ( REGION. EQ. 1 .OR.REGION.EQ. 3) GO TO 311 
NF=NF-1 
WF=WFSTD-DWF 
311 DO 60 1=1, NF 

COX=l. 

CFU=1 . 

L=1 

WTOT= WOX+WF 

OF=WOX/WF 

FRACT= 1 ./ (OF+1. ) 

CALL CHMBR (FRACT,CSTR,GC ,MW) 

PC= CSTR*WTOT/(AT*GCON) *CSTREF 
PCSI=PC/144 . 

B= (GC+1 . ) / (GC-1 . ) 

A= GC*SQRT ( (2 ./ (GC+1. ) ) **B) *CSTR 

FG= A*SQRT(SM**2.+(SL*PI*R/LL)**2.)/R 

FREQN=FG/(2.*PI) 

ALFAZ=0 . 

FREQZ=FG 

S (1) = CMPLX { 0 . ,FG) 



IF (INPUTF.EQ. 1) S (1) =CMPLX (FRINP , FIINP) 

IF (JP.GT.l.OR.I.GT.l) S(1)=SZ 
IF (M . GT . 1 . AND . I . EQ . 1 ) S ( 1 ) =S SAVE 
ALFG= REAL(S (1) ) 

FG=AIMAG(S (1) ) 

AA=A*A 

RHOC=GC * PC *GCO N/AA 
VZ=WTOT/(PI*RHOC*R**2. ) 

VOA=VZ/A 
W=VZ*VZ 
AMV=AA-W 
AMOR=AA* (SM/R) **2 . 

MACH=W/AA 

CALL WOP (GAMOX , CORFOX , WOX , AOX , TOMOX , VOX , LOX , COX , ROX , IOX , RHOX ,• 
*DELPO , IFLAGO , PCHOKO ) 

DELPOX=DELPO /144. 

PDOX= (PC+DELPO )/144. 

DPOPCO=DELPO /PC 

7 CALL WOP (GAMF ,CORFF , WF , AF ,TOMF , VF , LF , CFU , RFU , IFU , RHOF , DELPFU , 

* IFLAGF , PCHOKF ) 

DELPF=DELPFU/144 . 

PDF= (PC+DELPFU) /144 . 

DPOPCF=DELPFU/PC 

IF (MAP . EQ . 1 .AND . IM. EQ . 1) CALL ERMAP 
IM=0 

1 DF=DELX*FG*MACH 

DA=DELX* . 1*ALFG 
DAA=10.*DELX 
IF (DAA.GT.DA) DA=DAA 
DO 10 J=1 , 3 
NBOX=0 . 

NBF=0. 

IF ( IFLAGO. EQ.l) GO TO 4 

CALL NBS ( IOX , ROX , COX , TAUOX , S ( J ) , NBOX , VCTROX , THETAO ) 

4 IF (IFLAGF. EQ. 1) GO TO 8 

CALL NBS ( IFU , RFU ,CFU ,TAUF , S ( J) ,NBF , VCTRF ,THETAF) 

8. NB (J) = (WOX*NBOX+WF*NBF) /WTOT 

CALL NCS (S ( J) , NC ( J ) ) 

ERR ( J) =NB ( J) -NC ( J) 

IF (CABS (ERR(J) ) . LE . .01) GO TO 50 
ER ( J) =REAL (ERR ( J) ) 

El (J) =AIMAG (ERR (J) ) 

IF (J-2) 5,6,10 

5 FREQZ=FG+DF 

S ( 2 ) =CMPLX (ALFG , FREQZ ) 

GO TO 10 

6 ALFAZ=ALFG+DA 

S ( 3 ) =CMPLX ( ALF AZ , FG ) 

FREQZ=FG 
10 CONTINUE 

DERDF= (ER(2) -ER(1) ) /DF 
DERDA= (ER ( 3 ) -ER ( 1) ) /DA 
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DEIDF=(EI(2)-EI(1) )/DF 
DEIDA= (El (3) -El (1) ) /DA 

ALFAZ=ALFG- (El (1) /DEIDF-ER (1) /DERDF) / (DE I DA/DE IDF -DERDA/DERDF ) 
FREQZ=FG- (El (1) /DEIDA-ER (1) /DERDA) / (DE IDF /DEIDA-DERDF /DERDA) 
SZ=CMPLX (ALFAZ ,FREQZ) 

IF(L-30) 15,17,17 

17 WRITE (6,173) L,SZ,S(1),SS 

WRITE (6,170) ERR(J) , ESAVE2 , ESAVE1 
GO TO 55 

15 IF (ABS ( (ALFG-ALFAZ) /ALFAZ) .GT . . 02 .OR.ABS (ALFG-ALFAZ) 

*.GT.l.) GO TO 9 

IF (ABS ( (FG-FREQZ)/(FG*MACH) ) .LE. .05) GO TO 45 
9 FG=FREQZ 

ALFG= ALFAZ 
ESAVE2=ERR(J) 

IF (L.EQ. 28) ESAVEl=ERR (J) 

SS=S(1) 

S ( 1) = CMPLX( ALFAZ, FREQZ) 

L= L+l 
GO TO 1 

45 IF (IFLAGO.EQ.l) GO TO 46 

CALL NBS (IOX,ROX,COX,TAUOX,SZ ,NBOX,VCTROX,THETAO) 

46 IF (IFLAGF .EQ . 1) GO TO 47 

CALL NBS ( I FU , RFU , CFU , TAUF , S Z , NBF , VCTRF , THETAF ) 

47 CALL NCS (SZ ,NC ( J) ) 

GO TO 51 

50 SZ=S(J) 

51 CSROX=l./(COX*SZ) 

CSRF=1./(CFU*SZ) 

ISOX=IOX*SZ 
ISF=IFU*SZ 
TAUSOX=-TAUOX*SZ 
TAUSF=— TAUF*SZ 
TSROX=REAL (TAUSOX) 

TSRF=REAL (TAUSF) 

XTSROX=EXP (TSROX) 

XTSRF=EXP (TSRF) 

TS IOX=AIMAG (TAUSOX) 

TSIF=AIMAG (TAUSF) 

SZI=AIMAG(SZ) 

HZ=SZI/(2.*PI) 

WOWO=HZ/FREQN 

GROW=REAL(SZ) 

SZHZ=CMPLX (GROW , HZ ) 

AOW=GROW/HZ 

IF (L .GE . 30) SZ=SSAVE 

IF (I.EQ.l) SSAVE=SZ 

52 WRITE (6,159) 

WRITE (6,160) WOX , PDOX , DELPOX , ROX , COX , CSROX , DPOPCO 
WRITE ( 6 , 161) WF , PDF , DELPF , RFU , CFU , CSRF , DPOPCF 
WRITE (6,162) 

IF (IFLAGO .NE. 1) GO TO 63 
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WRITE (6,720) PCHOKO ,DPOPCO 
GO TO 64 

63 WRITE (6,168) IOX,ISOX,XTSROX,TSIOX,NBOX,VCTROX,THETAO 

64 IF (IFLAGF.NE. 1) GO TO 65 
WRITE (6, 7 21) PCHOKF ,CPOPCF 
GO TO 67 

65 WRITE (6,169) IFU , ISF ,XTSRF ,TSIF ,NBF ,VCTRF ,THETAF 

67 WRITE (6,163) 

WRITE (6,164) PCS I ,CSTR,A,VOA,OF ,NC (J) ,L 
WRITE (6,165) SZHZ ,FREQN , WOWO , AOW 
WRITE (6,166) 

55 WF=WF+DWF 

IF ( REGION. EQ . 2 .OR. REGION . EQ • 4) WF=WF-2 . *DWF 
60 CONTINUE 

WOX=WOX+DW 

IF (REGION .EQ . 1 .OR . REGION . EQ . 4 ) WOX=WOX-2 . *DW 
70 CONTINUE 

IF (WMAX .EQ. WMIN . AND .WFMAX . EQ . WFMIN) GO TO 500 
IF (WMAX . EQ . WOSTD . AND . REGION. EQ . 2) GO TO 500 
IF ( WFMAX. EQ. WFMIN. AND. REGION. EQ. 2) GO TO 500 
IF (REGION.LT. 4) GO TO 440 
500 CONTINUE 

1000 CONTINUE 

149 FORMAT ( 3F12 . 5 , 112 ,G12 . 5 , 112 ) 

150 FORMAT ( 6G12.5) 

151 FORMAT ( 7F10 . 5 ,F2 . 0) 

152 FORMAT ( 1H1 , 5X , 23HGAS-GAS STABILITY MODEL) 

153 FORMAT (1H0 , 27HINPUT PARAMETER FOR OXIDANT) 

154 FORMAT (1H ,1X,5HW MAX,7X,5HW MIN , 7X , 5HGAMMA, 7X, 7HORIF CF,5X, 
*9HORIF AREA, 3X , 6HLENGTH , 6X , 8HDOME VOL , 4X, 3HTAU , 9X, 4HT/MW, 

*9X, 2HDW, 9X , 5HW STD) 

155 FORMAT (1H ,11G12.5) 

156 FORMAT (1H0 , 24HINPUT PARAMETER FOR FUEL) 

157 FORMAT ( 1H0 , 18HCHAMBER PARAMETERS) 

158 FORMAT (1H , IX , 6HLENGTH , 6X , 6HRADIUS , 6X , 8HTHROAT A,4X,6HC* EFF , 

*6X, 6HWAVE M, 6X,6HWAVE L,6X,15HNOZZLE RESPONSE , 9X , 2HIP) 

159 FORMAT (1H0 , 12X , 1HW, 11X , 6HPD , PSI , 6X , 9HDEL P ,PSI , 3X ,1HR, 11X , 1HC , 

* 11X , 8HRE ( 1/CS ) , 4X, 8HIM (1/CS) , 4X , 7HDELP/PC) 

160 FORMAT (1H , 7HOXIDANT , 2X ,G12 . 4 ,F9 . 3 , 3X , F12 . 3 , 4G12 . 4 ,F12 . 7 ) 

161 FORMAT (1H , 4HFUEL , 5X,G12 . 4 ,F9 . 3 , 3X ,F12 . 3 , 4G12 . 4 ,F12 . 7) 

162 FORMAT (1H0 , 12X , 1HI , 11X, 6HRE (IS) ,6X,6HIM(IS) ,4X,11HE(RE (-T*S) ) , 
*2X, 10HIM (-TAU*S) , 3X , 8HRE (RESP) , 4X , 8HIM (RESP) , 4X , 9HVCTOR ,T=0 , 3X , 

* 9HTHETA , T=0 ) 

163 FORMAT (1H0 , 12X , 6HPC , PSI , 6X, 2HC* , 10X , 9HSOUND SPD, 3X, 4HMACH , 8X , 

* 3HO/F , 9X , 8HRE (RESP) , 4X , 8HIM (RESP) , 4X , 10HITERATIONS) 

164 FORMAT (1H , 7HCHAMBER, 2X,F9 . 3 , 3X , 7G12 . 4 ) 

165 FORMAT (1H0 , 3X, 12HGROWTH RATE= , F12 . 4 , 2X , 14HFREQUENCY (HZ ) = , 

* F10 . 3 , 2X , 22HNATURAL FREQUENCY (HZ )= , F10 . 3 , 2X, 5HW/WO= ,F8 . 5 , 

*2X, 4HA/W=,F10 .6) 

166 FORMAT (1H , 131H*************************************************** 
******************************************************************* 

***************) 
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168 FORMAT (1H , 7HOXIDANT , 2X, 9G12 . 4) 

169 FORMAT ( 1H , 4HFUEL , 5X, 9G12 .4) 

170 FORMAT (1H , 42HERROR IN RESPONSE FUNCTIONS ARE FOR L=30 , 2F9 . 4 , 
*6H L=29 , 2F9 . 4 , 9X , 6H L=28 ,2F9.4) 

171 FORMAT (1H0, 4 3HUNITS USED ARE LB-FT-SEC-DEG R UNLESS NOTED) 

172 FORMAT (1H0/1H0/1H0/1H0) 

173 FORMAT (1H , 30HABS (ERROR) DI VERGES. ITERATION=, 12 ,10H CALC' D S=, 
*2F9 . 2 , 6H S ( 1) = , 2F9 . 2 , 15H PREVIOUS S(1)=,2F9.2) 

174 FORMAT (4G12. 6, 312) 

175 FORMAT (2F 12. 6, 112) 

499 FORMAT (2F20. 10 ,212) 

708 FORMAT (2110) 

709 FORMAT (6E12.5) 

710 FORMAT (1H0, IX, 5HDEL X , 8X , 4HPMAX , 8X , 4HPMIN , 6X , 5HWHICH , 8X , 2HNP , 
*9X, 3HMAP) 

720 FORMAT (1H ,26HOXIDANT IS CHOKED, PCHOKE=,Fl2 .5,5X, 

*15HDELP CHOKED/PC= , F 1 2 . 5 ) 

721 FORMAT (1H ,23HFUEL IS CHOKED, PCHOKE= ,F12 . 5 , 8X , 

*15HDELP CHOKED/PC= , F 1 2 . 5 ) 

730 FORMAT (1H0 , IX , 21HINITIAL (GROWTH , FREQ) , 2X , 6HINPUTF) 

STOP 

END 

$ IBFTC WO 

SUBROUT INE WOP ( GC , C , W , A , T , V , D , AC , AR , AI , RHOD , DELP , IFLAG , PCHOKE ) 
COMMON / AMAIN/PI, GCON,GASC, PC, WTOT 
DIMENSION PDKNO (5) ,WKO(5) ,TERMl(5) ,TERM2(5) 

G=GC 
IFLAG=0 
GP= (G+l . ) /G 
GTW=2./G 
GM=(G-1.)/G 
GMR= l./GM 
GR=1./G 
TOM=T 

PCHOKE=PC/( (2. /(G+l.) ) **GMR) 

TERMO= (PC/PCHOKE) **GTW 
TERMT= (PC/PCHOKE) **GP 
PCGM=PC**GM 
PCGR=PC**GR 

WCHOKE=A*C*PCHOKE*SQRT ( 2 . *GCON*GMR* (TERMO-TERMT) / (GASC*TOM) ) 

IF( (W-WCHOKE) .GT.0.) GO TO 50 
DPD= ( PCHOKE-PC ) /4 . 

PDKNO (1)=PCH0KE 
PDKNO ( 2 ) =PCHOKE-DPD 
PDKNO (3) =PDKNO (2) -DPD 
PDKNO (4 )= PDKNO (3) -DPD 
PDKNO ( 5) = PC 
WKO ( 1 ) =WCHOKE 
WKO (5) =0 . 

DO 10 1=2,4 

TERMl (I ) = ( PC/PDKNO (I) ) * *GTW 
TERM2 (I) = (PC/PDKNO (I) ) **GP 
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WKO (I) =A*C*PDKNO (I) *SQRT (2 . *GCON*GMR* (TERM1 (I) -TERM2 (I) ) / 

* (GASC*TOM) ) 

10 CONTINUE 

ONE= (W-WKO ( 2 ) ) * (W-WKO (3) ) * (W-WKO (4) ) * (W-WKO(5) )/ 

* ( (WKO (1) -WKO (2) ) * (WKO ( 1) -WKO (3) ) * (WKO(l) -WKO (4) ) * (WKO(l) -WKO (5) ) ) 
TWO= (W-WKO (1) ) * (W-WKO (3) ) * (W-WKO (4) ) * (W-WKO (5) )/ 

* ( (WKO (2) -WKO(l) )* (WKO (2) -WKO (3) ) * (WKO (2) -WKO (4) ) * (WKO (2) -WKO (5) ) ) 
THR= (W-WKO (1) ) * (W-WKO (2) ) * (W-WKO (4) ) * (W-WKO (5) )/ 

* ( (WKO (3) -WKO (1) ) * (WKO (3) -WKO (2) ) * (WKO (3) -WKO (4) ) * (WKO (3) -WKO (5) ) ) 
FOR= (W-WKO (1) ) * (W-WKO (2) ) * (W-WKO (3) ) * (W-WKO (5) )/ 

* ( (WKO (4) -WKO(l) )* (WKO (4) -WKO (2) ) * (WKO (4) -WKO (3) ) * (WKO (4) -WKO (5) ) ) 
FI V= (W-WKO ( 1 ) ) * (W-WKO ( 2 ) ) * (W-WKO ( 3 ) ) * (W-WKO (4 ) ) / 

* ( (WKO (5) -WKO(l) )* (WKO (5) -WKO (2) ) * (WKO (5 ) -WKO (3) ) * (WKO (5) -WKO (4) ) ) 
PD=PDKNO (1) *ONE+PDKNO (2) *TWO+PDKNO ( 3) *THR+PDKNO (4) *FOR+PDKNO (5) * 
*FIV 

20 RHOD=PD/(GASC*TOM) 

DELP=PD-PC 
GMM=(1.-G)/G 
PCGMM=PC**GMM 
PDGMM=PD * * GMM 
PDGM=PD* *GM 

ALPHA=1 . /G+GM/2 . * 1 . / ( 1 . -PCGMM*PDGM) 

BETA=GM/2 . * 1 . / ( 1 . -PCGM*PDGMM) 

AR=-1. /ALPHA 
AI=W*D/ (GCON*PC*A) 

AC=RHOD*V/ (BETA*G*W) 

AC =- AC* ALPHA 
GO TO 60 

50 TERMTH=SQRT (GMR* (TERMO-TERMT ) *2 . *GCON/ (GASC*TOM) ) 

PD=W/ (A*C*TERMTH) 

IFLAG=1 
DELP=PD-PC 
60 RETURN 

END 
$IBFTC NB 

SUBROUTINE NBS ( AI , AR , AC , TAU , S , NB , VECT , THET ) 

COMPLEX S , EXPO , NB , NOTAU 
EXPO=CEXP (-TAU*S) 

NOTAU=-AC*S / (1 .+AC*AR*S+AC*AI*S*S) 

NB= NOTAU*EXPO 
VECT=CABS (NOTAU) 

THET=ATAN2 (AIMAG (NOTAU) , REAL (NOTAU) ) *180 ./ 3 .1416 

RETURN 

END 

$IBFTC CH 

SUBROUTINE CHMBR (XA ,CA ,GA , A) 

COMMON /BF/F 

DIMENSION F ( 200) ,B(50,4) ,T(3) ,U(3,3) ,V(3) 

EQUIVALENCE (F ( 1) ,B(1,1) ) 

X= XA 
1=0 

10 1= 1+1 
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IF (B(I,1) . LT . X ) GO TO 10 
T ( 2) = B ( 1 , 1) 

U ( 2 , 1) = B (I , 2) 

U ( 2 , 2) = B (1 , 3) 

U (2 , 3) = B (1 , 4 ) 

1= 1+1 
T(3)= B ( 1 , 1 ) 

U (3 , 1) = B (1 ,2) 

U (3 , 2) = B(I,3) 

U ( 3 , 3) = B(I,4) 

1= 1-2 

T (1) = B (1 , 1) 

U(l,l)= B (1 , 2) 

U (1 , 2) = B (1 , 3) 

U (1 , 3) = B(I,4) 

ZERO= (X-T ( 2) ) * (X-T ( 3) ) / ( (T (1) -T (2) ) * (T (1) -T ( 3 ) ) ) 
ONE = (X-T (1) ) * (X-T (3))/((T(2) -T (1) ) *(T(2) -T (3) ) ) 
TWO = {X-T (1) ) * (X-T (2) )/((T(3)-T(l) ) * (T (3) -T (2) ) ) 
DO 20 M=1 , 3 

V (M) = U (1 ,M) *ZERO+U ( 2 ,M) *ONE+U (3 ,M) *TWO 
20 CONTINUE 
CA= V(l) 

GA= V ( 2) 

A= V (3) 

RETURN 
END 
$IBFTC NS 

SUBROUTINE NCS(S,NC) 

COMPLEX EXPT , EXPO , NC , BRAKO , BRAKT 

COMPLEX S , GN , SQR , XONE , XTWO , D , C , E , XONEL , XTWOL 

REAL LL 

COMMON /SCLL/VZ ,AMOR,AMV,GC ,W ,LL,AA,GNR,GNI 
GN=CMPLX ( GNR , GN I ) 

SQR=CSQRT( (VZ*S*VZ*S)+ (S*S+AMOR) *AMV) 

XONE= (VZ*S+SQR)/AMV 
XTWO= (VZ*S-SQR)/AMV 
D= AMV+GC*GN*W 
E= (GC*GN-1.)*VZ*S 
C=- (XONE*D+E) / (XTWO*D+E) 

XONEL=-XONE*LL 
XTWOL=- XTWO * LL 
EXPO=CEXP (XONEL) 

EXPT=CEXP (XTWOL) 

BRAKO=XONE*EXPO+XTWO*C*EXPT 

BRAKT=EXPO+C*EXPT 

NC=1 ./GC-AA* BRAKO/ (GC*VZ* (S*BRAKT+VZ*BRAKO) ) 

RETURN 

END 

$IBFTC BD 

BLOCK DATA 
DIMENSION F (200) 

COMMON /BF/F 


COMMON /AMAIN/PI ,GCON,GASC, PC, WTOT 
DATA PI ,GCON ,GASC/3 .1416,32.17,1546./ 

DATA (F ( J) , J=1 , 50 ) / . 004 , .006, .008, .01042, .010989, .011628, .012346 

A. 013158, .014085, .015152, .016 39 6, .017857, .019608, .021739 , .024390, 

B. 027778, .032258, . 038462 ,. 047619 , .062500, . 0666667 071429 , .076923, 

C. 08 3333, .090909, . 111111 ,. 14 2857 , .200000, .208333, .217391, .227273, 

D. 238095, .250000, . 263158 ,. 277778 ,. 294118 , .312500, .333333, .350000, 

E. 400000 , .450000, .500000, . 550000,-60, . 65 , . 70 , . 75 , . 80 , . 85 , .90/ 

DATA (F(J) ,J=51,100)/2210. , 

A2516. ,2781. ,3062. ,3124. ,3194. , 3265 ., 3345 ., 3433 ., 3532 . ,3641. ,3764. 
B3904. ,4064. ,4249. ,4462. ,4712. ,5008 . ,5374 . ,5867. ,5991. ,6127. ,6277. 
C6445 . ,6634. ,7089. ,7653. ,8216. ,8261. ,8300. ,8335. ,8363. ,8384. ,8397. 
D8402 . ,8398. ,8384. ,8358. ,8332. ,8222. ,8077. ,7903. ,7705. ,7491. ,7263. 
E7023. ,6773. ,6511. ,6237 ., 59 47 ./ 

DATA (F ( J) ,J=101,150)/1. 3259, 1.3080, 1.2958, 1.2829, 1.2801, 

Al. 2769, 1.2734, 1.2694, 1.264 8, 1.259 4, 1.2529, 1.2449, 1.2350, 1.2227, 
Bl. 2078, 1.1907, 1.1727, 1.1559 ,1.1421,1.1321,1.1306,1.1292,1.1280, 
Cl. 1270, 1.1261, 1.1254, 1.13 31, 1.1650, 1.1716, 1.1791, 1.1876, 1.1971, 
Dl. 2077, 1.2191, 1.2312, 1.2436 ,1.2559 ,1.2681,1.2766,1.2986,1.3177, 
El. 3353, 1.3517, 1.3661, 1.376 3, 1.3842, 1.3893, 1.3922, 1.3943, 1.3966/ 
DATA (F(J) ,J=151,200)/ 

A31. 136, 30. 721, 30. 318, 29. 844, 29. 733, 29. 612, 29. 476, 29. 323, 29. 151, 
B28. 956, 28. 731, 28. 469, 28. 160, 27. 789, 27. 732, 26. 751, 25. 984, 24 .931, 
C23. 430, 21. 210, 20. 642, 20. 020, 19. 337, 18. 585, 17. 750, 15. 772, 13. 225, 
D9. 946, 9. 578, 9. 20 3, 8. 823, 8. 4 37, 8. 046, 7. 651, 7. 252, 6. 852, 6. 450, 6. 048 
E5. 760, 5. 040, 4. 480, 4. 032, 3. 665, 3. 360, 3. 102, 2. 880, 2. 688, 2. 520, 2. 372 
F2.240/ 

END 
$IBFTC RS 

SUBROUTINE RSTRP(I,J) 

REAL LL , LOX , LF 
DIMENSION V ( 14 ) 

COMMON/ERM/FRMAX ,FIMAX ,FRMIN ,FIMIN ,FRCUT ,FICUT , IOX ,ROX ,COX , 

*TAUOX , CFU , RFU , IFU , TAUF , IFLAGO , IFLAGF , WF , WOX 
COMMON/WH/ R, AT ,TOMOX,TOMF , AOX ,AF ,LOX ,LF ,VOX ,VF 

COMMON /SCLL/VZ , AMOR , AMV , GC , W,LL,AA,GNR,GNI 

IF(I.GT.l) GO TO 5 
V (1) =LL 
V(2)=R 
V ( 3 ) =AT 
V ( 4 ) =TOMOX 
V ( 5 ) =TOMF 
V (6) =AOX 
V ( 7 ) =AF 
V ( 8 ) =LOX 
V ( 9 ) =LF 
V ( 10 ) =VOX 
V ( 1 1 ) =VF 
V (12) =TAUOX 
V ( 1 3 ) =TAUF 
RETURN 

GO TO (11,12,13,14,15,16,17,18,19,20,21,22,23) ,J 
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11 LL=V (1 ) 

RETURN 

12 R=V ( 2 ) 

RETURN 

13 AT=V ( 3 ) 

RETURN 

14 TOMOX=V ( 4 ) 

RETURN 

15 TOMF=V ( 5 ) 

RETURN 

16 A0X=V(6) 

RETURN 

17 AF=V ( 7 ) 

RETURN 

18 LOX=V ( 8 ) 

RETURN 

19 LF=V(9) 

RETURN 

20 VOX=V(10) 

RETURN 

21 VF=V (11) 

RETURN 

22 TAUOX=V (12) 

RETURN 

23 TAUF=V (13) 

RETURN 

END' 

$IBFTC WHH 

SUBROUTINE WHICHP (WHICH ,JP ,PMAX , DP) 

REAL LL , LOX , LF 
INTEGER WHICH 

COMMON /SCLL/VZ , AMOR,AMV ,GC , W ,LL ,AA ,GNR ,GNI 

COMMON/ERM/FRMAX ,FIMAX ,FRMIN , FIMIN , FRCUT ,FICUT , IOX , ROX , COX , 
*TAUOX , CFU , RFU , IFU , TAUF , IFLAGO , I FLAG F , WF , WOX 
COMMON/WH/ R , AT , TOMOX , TOMF , AOX ,AF ,LOX ,LF ,VOX ,VF 

GO TO (410,411,412,413,414,415,416,417,418,419,420,421, 
*422) , WHICH 

410 IF(JP.EQ.l) LL=PMAX+DP 
LL=LL-DP 

WRITE (6,480) LL 
GO TO 400 

411 IF ( JP .EQ . 1) R =PMAX+DP 
R=R-DP 

WRITE (6 , 481) R 
GO TO 400 

412 IF (JP.EQ. 1)AT =PMAX+DP 
AT=AT-DP 

WRITE (6,482)AT 
GO TO 400 

413 IF (JP.EQ, l)TOMOX=PMAX+DP 
TOMOX=TOMOX-DP 

WRITE (6, 483) TOMOX 



414 


GO TO 400 

IF ( JP .EQ . 1) TOMF =PMAX+DP 
TOMF=TOMF-DP 
WRITE (6, 484) TOMF 
GO TO 400 

415 IF { JP . EQ . 1 ) AOX =PMAX+DP 
AOX=AOX-DP 

WRITE (6,485) AOX 
GO TO 400 

416 IF ( JP .EQ . 1) AF =PMAX+DP 
AF=AF-DP 


WRITE (6, 486) AF 
GO TO 400 

417 IF (JP.EQ.l)LOX =PMAX+DP 

LOX=LOX-DP 

WRITE (6, 487) LOX 
GO TO 400 

418 IF ( JP . EQ .1) LF =PMAX+DP 

LF=LF-DP 

WRITE (6,488) LF 
GO TO 400 

419 IF (JP.EQ.l) VOX =PMAX+DP 

VOX=VOX-DP 

WRITE (6,489) VOX 
GO TO 400 

420 IF (JP.EQ. 1)VF =PMAX+DP 

VF=VF-DP 


WRITE (6,490) VF 
GO TO 400 

421 IF (JP.EQ. l)TAUOX=PMAX+DP 
TAUOX=TAUOX-DP 
WRITE (6,491) TAUOX 
GO TO 400 

ITAUF =PMAX+DP 


422 

IF ('JP . EQ . 1 
TAUF=TAUF- 
WRITE (6,49 

480 

FORMAT ( 1H 

481 

FORMAT (1H 

482 

FORMAT (1H 

483 

FORMAT (1H 

484 

FORMAT (1H 

485 

FORMAT (1H 

486 

FORMAT (1H 

487 

FORMAT ( 1H 

488 

FORMAT (1H 

489 

FORMAT (1H 

490 

FORMAT ( 1H 

491 

FORMAT (1H 

492 

FORMAT (1H 

400 

RETURN 

END 

$ IBFTC MP 


, 15HCHAMBER LENGTH= , 

, 15HCHAMBER RADIUS=, 

, 19HNOZZLE THROAT AREA= , 

, 20HOXIDANT TEMP/MOL WT=, 

, 20HFUEL TEMP /MOL WT= , 

, 21HOXIDANT ORIFICE AREA= , 

, 18HFUEL ORIFICE AREA= , 

, 24HOXIDANT INJECTOR LENGTH= , 

, 21HFUEL INJECTOR LENGTH= , 

, 20HOXIDANT DOME VOLUME=, 

, 17HFUEL DOME VOLUME=, 

, 30HOXIDANT COMBUSTION DELAY TIME=, 
, 27HFUEL COMBUSTION DELAY TIME= , 


5X ,G11 . 5 ) 
5X,Gll .5) 
5X ,G11 . 5) 
5X ,G11 . 5 ) 
5X ,G11 . 5) 
5X ,G11 . 5) 
5X ,G11 . 5 ) 
5X ,G11 . 5) 
5X ,G11 . 5) 
5X,G11 . 5) 
5X,G11. 5) 
5X ,G11 . 5) 
5X ,G11 . 5) 


SUBROUTINE ERMAP 
REAL LL , IOX , IFU 

INTEGER FRCUT , FICUT , QUAD (20,40) 

COMPLEX S, NB , NC , NBOX , NBF , DRESP 
DIMENSION ERR (20,40) ,TR(20) ,TI (40) 

COMMON /SCLL/VZ ,AMOR, AMV ,GC , W,LL,AA,GNR,GNI 

COMMON/E RM/FRMAX , F IMAX , FRMI N , F IMIN , FRCUT , F I CUT , IOX , ROX , COX , 
*TAUOX , CFU , RFU , IFU , TAUF , IFLAGO , IFLAGF , WF , WOX 
COMMON /AMAIN/PI, GCON,GASC, PC, WTOT 
FR=FRMAX 
DFR=0 . 

DFI=0. 

IF (FRCUT . NE . 1) DFR= (FRMAX-FRMIN) /FLOAT (FRCUT) 

IF (FICUT .NE .1) DFI= (FIMAX-FIMIN) /FLOAT (FICUT) 

DO 1000 1=1, FRCUT 
FI=FIMAX 

DO 1010 J=l, FICUT 
S=CMPLX (FR, FI ) 

IF ( IFLAGO. EQ.l) GO TO 40 

CALL NBS (IOX,ROX,COX,TAUOX,S ,NBOX,A,B) 

40 IF (IFLAGF. EQ.l) GO TO 80 

CALL NBS ( IFU , RFU , CFU , TAUF , S ,NBF , A ,B) 

80 NB=(WOX*NBOX+WF*NBF)/WTOT 

CALL NCS (S ,NC) 

DRESP=NB-NC 

ANGLE=ATAN2 (AIMAG (DRESP) , REAL (DRESP) ) 

IF (ANGLE) 500,510,520 
500 IF (ANGLE+PI/2 . ) 600,601,601 
510 QUAD (I , J) =0 

GO TO 90 

520 IF (ANGLE-PI/2.) 605,605,606 

600 QUAD (I, J) =3 
GO TO 90 

601 QUAD (I, J) =4 
GO TO 90 

605 QUAD ( I , J ) = 1 
GO TO 90 

606 QUAD (I, J) =2 
GO TO 90 

90 ERR ( I , J) =CABS (NB-NC) 

FI=FI-DFI 
1010 CONTINUE 
TR (I ) =FR 
FR=FR-DFR 
1000 CONTINUE 
FI=FIMAX 

DO 1020 J=l, FICUT 
TI (J) =FI 
FI=FI-DFI 
1020 CONTINUE 

WRITE (6, 9 20) 

WRITE (6,900) (TR ( J) ,J=1, FRCUT) 



WRITE (6,905) 

DO 1030 J=1,FICUT 

WRITE (6,910) TI(J) , (ERR(I,J) ,QUAD(I,J) ,I=1,FRCUT) 
1030 CONTINUE 

WRITE (6,905) 

900 FORMAT ( 1H , IX , 4HFREQ , 58X, 6HGROWTH/1H ,8X,15F8.1) 

905 FORMAT (1H ,128H 

•k _ _ . . . 


* ) 

910 FORMAT (1H ,F6 . 0 , 1H- , IX , 15 (F6 . 2 , 1H* , II) ) 

920 FORMAT (1H0, 5 8HMAP OF ERRORS, WHICH IS ABS (NB-NC) 

* ERROR) 

RETURN 

END 


$DATA 

2.658 2.658 

. 36574E-02 .50255E-04 

.505 .3 

. 24363E-02 .50255E-04 

.5 .1575 .0185 

2.658 .505 

. 19514E-03 .001 

750. -750. 90000 


1.36 

1.41 


.1 


16.875 

393. 

1 . 


-5. 

50000 


1. . 19514E-02 

1 

1. . 30957E-02 

1 

.9166 0. 

.2 

9900. 621 

1540 


AND * QUADRANT OF 


. 27833E-03 
.5 

. 98092E-02 
1.840. 
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APPENDIX B 


PROGRAM INPUT FORMAT 


Card 1 6G12.5: 

1-12 WMAX 
13-24 W3VDN 

25-36 GAMOX 
37-48 CORFOX 
49-60 AOX 
61-72 LOX 
Card 2 3G12.5, 112, 
1-12 VOX 
13-24 TAUOX 
25-36 TOMOX 
37-48 IP 

49-60 DELX 
Card 3 6G12. 5: 

1-12 WFMAX 
13-24 WFMIN 

25-36 GAMF 
37-48 CORFF 
49-60 AF 
61-72 LF 

Card 4 3G12.5, 112: 
1-12 VF 
13-24 TAUF 
25-36 TOMF 
37-48 MAP 


maximum oxidant flow rate, lb/sec 

minimum oxidant flow rate, lb/sec 

(= WMAX if no variation in oxidant flow rate) 

oxidant specific heat ratio 

oxidant injector orifice coefficient 

o 

total oxidant injector area, ft 
oxidant injector length, ft 
G12.5: 

3 

total oxidant dome volume, ft 

oxidant combustion delay time , sec 

oxidant temperature/molecular weight, °R-mole/lb 

number of cases in parametric study 

(= 0 or 1 if only one case) 

constant used to increment S (e. g. , DELX = 0. 5) 

maximum fuel flow rate , lb /sec 
minimum fuel flow rate, lb/sec 
(= WFMAX if no variation in fuel flow rate) 
fuel specific heat ratio 
fuel injector orifice coefficient 

2 

total fuel injector orifice area, ft 
fuel injector length, ft 

3 

total fuel dome volume, ft 

fuel combustion delay time, sec 

fuel temperature/molecular weight, °R -mole/lb 

= 1 if error map desired 

= 0 if no error map 
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Card 5 7F10. 5, F2.0: 

1-10 LL combustion chamber length, ft 

11-20 R chamber radius, ft 

o 

21-30 AT nozzle throat area, ft 

31-40 CSTREF C* efficiency 
41-50 GNR real part of nozzle admittance 
51-60 GNI imaginary part of nozzle admittance 

61-70 SM transverse wave mode 

71-72 SL longitudinal wave mode 

Card 6 4G12.6: 

1-12 WOSTD standard oxidant flow rate, lb/sec 

(WOSTD = WMAX if no variation in oxidant rate) 
13-24 WFSTD standard fuel flow rate, lb/sec 

(WFSTD = WFMAX if no variation in fuel rate) 
25-36 DW increment in oxidant rate , lb/sec 

37-48 DWF increment in fuel rate, lb/sec 
Card 7 (If IP = 0) 2G12. 6, 112: 

1-12 FRINP I initial guess of frequency (FRINP + iFIINP); 

13-24 FIINP J dimensions are 1/sec and Hz, respectively 
25-36 INPUT F if = 0 ignore above and uses natural frequency 
if = 1 use above guessed frequency 
Card 7 and on (If IP * 0) 4G12. 6, 312, Number of cards I = IP: 

1-12 Z1(I) parameter maximum value 

13-24 Z2(I) parameter minimum value 

25-36 Z3(I) 1 equivalent to FRINP and FIINP; 

37-48 Z4 (I) J units are 1/sec and Hz, respectively 

49-50 IZ1(I) = 0 to 14 (see note (1)) 

51-52 IZ2(I) number of increments between Z1(I) and Z2(I) 

53-54 IZ3(I) equivalent to INPUT F 

Last card (If MAP - 1) 4G12. 6, 212: 

1-12 FRMAX 

13-24 FRMIN't specify range of complex frequency values 
25-36 FIMAX J for error map; units are 1/sec and rad/sec 
37-48 FIMIN 

49-50 FRCUT number of increments between FRMAX and FRMIN 
with maximum = 15 

51-52 FICUT number of increments between FIMAX and FIMIN 
with maximum = 40 
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Note (1) - Parameter codes IZ1(I) are assigned as follows: 
IZ1(I) Parameter 


1 

LL 

chamber length 

2 

R 

chamber radius 

3 

AT 

nozzle area 

4 

TOMOX 

oxidant temperature/molecular weight 

5 

TOMF 

fuel temperature/molecular weight 

6 

AOX 

oxidant injector area 

7 

AF 

fuel injector area 

8 

LOX 

oxidant injector length 

9 

LF 

fuel injector length 

10 

VOX 

oxidant dome volume 

11 

VF 

fuel dome volume 

12 

TAUDX 

oxidant combustion delay time 

13 

TAUF 

fuel combustion delay time 
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Card 


Columns 



1 to 12 

13 to 24 

25 to 36 

37 to 48 

49 to 60 

61 to 72 

III 

2.658 

2.658 

1.36 

1 . 

0. 19514E-02 

0. 27833E-03 

1 - 

. 36574E-02 

. 50255E-04 

16.875 

1 

.5 



.505 

.3 

1.41 

1 . 

. 30957E-02 

.98092E-02 

Kg 

. 24363E-02 

. 50255E-04 

393. 

1 




Card 


Columns 



1 to 10 

11 to 20 

21 to 30 

31 to 40 

41 to 50 

51 to 60 

61 to 
70 

71 
& 

72 

5 

0.5 

0.1575 

0.0185 

1 . 

0.9166 

0.0 

1.84 

0 . 


Card 


Columns 


1 to 12 


13 to 24 


25 to 36 


37 to 48 
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fc 
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Figure 4. - Concluded. 
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